#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>
#include <vector>
#include <fstream>
#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include "math_functions.h"
#define windowSize 500
#define threadsPerBlock 512
#define EDGE 5
#define ae_pi 3.1415926535897932384626433832795
#define minDIV 0.0001
#define numDIV  600
#define FACTOR 10000
#define sizeAlpha 80
using namespace std;

double locallngamma(double x, double* sgngam);
void localAEPD_Fit(double * x, int N, double * params);
void MakeDeviationMatrix(double **D);

/////////////////////////////////////////////////////////////////////////////
// Code imported from CPUTime.cpp:
// This code is used to check performance
//////////////////////////////////////////////////////////////////////////////

typedef struct timeval TimeVal, *TimeValPtr;
/* timeval is a struct with fields tv_sec for seconds and
                                   tv_usec for microseconds */

int GetTime (TimeVal &pTime)
{
  struct rusage usage;
  int status;

  if ((status= getrusage(RUSAGE_SELF,&usage))!=0)
    return (status);  /* return error code */

  pTime=usage.ru_utime; /* user time */

  return (status);
}

long GetTimeDifference (TimeVal &pStart, TimeVal &pStop)
{
  long usec, sec;
  std ::cout << "In Get Time Difference "<< std::endl;
  sec=pStop.tv_sec-pStart.tv_sec;    /* difference in seconds */
  usec=pStop.tv_usec-pStart.tv_usec; /* difference in microseconds */


  return ( sec*1000000l +  usec); /* result in microseconds */
}

////////////////////////////////////////////////////////////////////////
// End of code to check time
////////////////////////////////////////////////////////////////////////

int main(void){

	  TimeVal startTime, endTime;
	  long totalTime;
	  ///////////////////////////////////////////////////////////////////////////////
	  // alpha and log-gamma(1/alpha) values are constant, declare and copy them here
	  ///////////////////////////////////////////////////////////////////////////////
	  int lengthDevAlpha = 80;
	  double *alpha, *logGammaofAlpha;
	  alpha = (double *) malloc (lengthDevAlpha*sizeof(double));
	  logGammaofAlpha = (double *) malloc (lengthDevAlpha*sizeof(double));

	  double sgngam = 0;

	  for (int i = 0; i < lengthDevAlpha; i++){
		  alpha[i] = 0.7 + i*0.1;
		  logGammaofAlpha[i] = locallngamma(1/alpha[i], &sgngam);
	  }
	  double *dev_alpha, *dev_lnGammaofAlpha;

	  cudaMalloc( (void **)&dev_alpha, sizeof(double) * lengthDevAlpha );
	  cudaMalloc( (void **)&dev_lnGammaofAlpha, sizeof(double) * lengthDevAlpha );

	  cudaMemcpy( dev_alpha, alpha, sizeof(double) * lengthDevAlpha, cudaMemcpyHostToDevice );
	  cudaMemcpy( dev_lnGammaofAlpha, logGammaofAlpha, sizeof(double)*lengthDevAlpha, cudaMemcpyHostToDevice );

	 //////////////////////////////////////////////////////////////////////
	 // Test Aleks' code with one launch
	 // First read in the data array (500 points) and then pass
	 //
	 /////////////////////////////////////////////////////////////////////

		double *dataArray, *param;
		dataArray = (double *) malloc (windowSize*sizeof(double));
		param = (double *) malloc (4 * sizeof(double));
		fstream myfile("gbpusd.txt", ios::in);
		if (!myfile.is_open()){
			cout << "Error opening gbpusd text file " << endl;
			return 1;
		}
		double d = 0;
		int index= 0;
		while (index < windowSize){
			myfile >> d;
			dataArray[index++] = d;
		}
		// Finished reading and assigning data
		/////////////////////////////////////////////////////////////////////////
		// Call C++ Version of AEPD_fit()
		/////////////////////////////////////////////////////////////////////////

		localAEPD_Fit(dataArray, windowSize, param);

		cout<< "Parameters from localAEPD_Fit" << endl;
		for (int i = 0; i<4; i++){
			cout << "Param[" << i << "] = " << param[i]<<endl;
		}

		////////////////////////////////////////////////////////////////////////////
		// Generate a deviation matrix:
		////////////////////////////////////////////////////////////////////////////
//		double deviationMatrix[numDIV][sizeAlpha];
//		double inc_A = 0.1;//////Iteration Increment of Alpha Parameter
//		double start_A = 0.7;//////// Value to Start Alpha Iteration
//		double *alpha, *deviationMatrix, *DIV;
//		alpha = (double *) malloc (sizeAlpha*sizeof(double));
//		deviationMatrix = (double *) malloc(sizeAlpha*numDIV*sizeof(double));
//		DIV = (double *) malloc(numDIV * sizeof(double));
/*		for (int i = 0; i < sizeAlpha; i++){
			  alpha[i] = start_A + i*inc_A;
		}

		for(int i = 0; i<numDIV; i++)
		{
			DIV[i] = minDIV*i;
		}

		for(int j = 0; j<numDIV; j++)
		{
			for(int i = 0; i<sizeAlpha; i++)
			{
				deviationMatrix[i + j*sizeAlpha] = pow(DIV[j],alpha[i]);

			///	fprintf(dMat,"%lf\t",D[j][i]);
			}
			///fprintf(dMat,"\n");
		}
		///fclose(dMat);


		///////////////////////////////////////////////////////////////////////////
		// Debugging code for deviation matrix: debug for one value of DIV = 0.0001
		///////////////////////////////////////////////////////////////////////////

		double tmpDIV;
		for (int i = 0; i< sizeAlpha; i++){
			tmpDIV = pow(minDIV, 0.7 + 0.1 * i);
			cout << "D[1][" << i << "] = "<< deviationMatrix[1*sizeAlpha + i] << "      tmpDIV = " << tmpDIV << endl;
		}
*/
	    myfile.close();
		return 0;
}





/*************************************************************************************
 * Log- Gamma function adapted from Alglib project
 *
 * Copy right and comments from original:
 ***************************************************************************************/

double locallngamma(double x, double* sgngam)
{

    double a;
    double b;
    double c;
    double p;
    double q;
    double u;
    double w;
    double z;
    int i;
    double logpi;
    double ls2pi;
    double tmp;
    double result;

    *sgngam = 0;

    *sgngam = 1;
    logpi = 1.14472988584940017414;
    ls2pi = 0.91893853320467274178;
    if( x<-34.0 )
    {
        q = -x;
        w = locallngamma(q, &tmp);
        p = floor(q);
        i = floor(p+0.5);
        if( i%2==0 )
        {
            *sgngam = -1;
        }
        else
        {
            *sgngam = 1;
        }
        z = q-p;
        if( z>0.5 )
        {
            p = p+1;
            z = p-q;
        }
        z = q*sin(ae_pi*z);
        result = logpi-log(z)-w;
        return result;
    }
    if( x<13 )
    {
        z = 1;
        p = 0;
        u = x;
        while(u>=3)
        {
            p = p-1;
            u = x+p;
            z = z*u;
        }
        while(u<2)
        {
            z = z/u;
            p = p+1;
            u = x+p;
        }
        if( z<0 )
        {
            *sgngam = -1;
            z = -z;
        }
        else
        {
            *sgngam = 1;
        }
        if( u==2)
        {
            result = log(z);
            return result;
        }
        p = p-2;
        x = x+p;
        b = -1378.25152569120859100;
        b = -38801.6315134637840924+x*b;
        b = -331612.992738871184744+x*b;
        b = -1162370.97492762307383+x*b;
        b = -1721737.00820839662146+x*b;
        b = -853555.664245765465627+x*b;
        c = 1;
        c = -351.815701436523470549+x*c;
        c = -17064.2106651881159223+x*c;
        c = -220528.590553854454839+x*c;
        c = -1139334.44367982507207+x*c;
        c = -2532523.07177582951285+x*c;
        c = -2018891.41433532773231+x*c;
        p = x*b/c;
        result = log(z)+p;
        return result;
    }
    q = (x-0.5)*log(x)-x+ls2pi;
    if(x>100000000 )
    {
        result = q;
        return result;
    }
    p = 1/(x*x);
    if( x>=1000.0 )
    {
        q = q+((7.9365079365079365079365*0.0001*p-2.7777777777777777777778*0.001)*p+0.0833333333333333333333)/x;
    }
    else
    {
        a = 8.11614167470508450300*0.0001;
        a = -5.95061904284301438324*0.0001+p*a;
        a = 7.93650340457716943945*0.0001+p*a;
        a = -2.77777777730099687205*0.001+p*a;
        a = 8.33333333333331927722*0.01+p*a;
        q = q+a/x;
    }
    result = q;
    return result;

} // end locallngamma

//////////////////////////////////////////////////////////////////////////////////////////////////////
//
// localAEPD_fit
// This function mimics AEPD_Fit(), has no external calls, this function will iterate over all alpha values
//
//////////////////////////////////////////////////////////////////////////////////////////////////////

void localAEPD_Fit(double * x, int N, double * params){

	// Take the array in x[0] - x[N-1] and put that into the vector data.
	vector<double> data(x, x+N);

	/*for (int i = 0; i < N; i++){
//		cout << "X[" << i << "] = " << x[i]<< endl;
		cout << "Before Sort data[" << i << "] = "<< data[i]<<endl;
	}*/


	/* Sorts the data: can be done in InArray.updateArray */
	double tmp;
	for(int i = 0;i < N-1 ;i++){
		for (int j = i + 1;j < N;j++){
			if (data[i] > data[j]){
				tmp  = data[i];
				data[i] = data[j];
				data[j] = tmp;
			}
		}
	}


	///sort(data.begin(),data.end());

	//for (int i = 0;i< N; i++){
	//	cout << "data sorted element " << i << " = " << data[i] << endl;
	//}

	//for (int i = 0;i< 4; i++){
	//	cout << "params " << i << " = " << params[i] << endl;
	//}
	double INCREMENT_ALPHA = 0.1;
	int sizeALPHA = 80;
	int NUM_STEPS = 500;

	double inc_A = INCREMENT_ALPHA;//////Iteration Increment of Alpha Parameter
	double start_A = 0.7;//////// Value to Start Alpha Iteration
	int sizeA = sizeALPHA;//// Number of iterations to perform

	double ** Kappa = new double*[sizeA];
	double ** Sigma = new double*[sizeA];
	double ** H = new double*[sizeA];
	for (int i = 0 ; i< sizeA ; i++){
		Kappa[i] = new double[N];
		Sigma[i] = new double[N];
		H[i] = new double[N];
	}

	double * X_minus = new double[sizeA];
	double * X_plus = new double[sizeA];


	int idx_k = 0;
	int idx_q = 0;
	double min_H = 10000000000;

	//int stepX = round(N/NUM_STEPS);
	int stepX = N/NUM_STEPS;
	int X_right = 0;
	int X_left = 0;
	if (stepX <1)
		stepX=1;

	for(int k=EDGE-1 ; k <= N-EDGE-1; k=k+stepX)
	{

		double p1 = data[k];

		X_right = N- k;
		X_left = k ;

		double * x_minus_tmp = new double[X_left];
		double * x_plus_tmp = new double[X_right];

		/* Remember that data[] is already sorted.
		 * Take k as the current index (p1 = data[k]),
		 * x_minus_tmp[0 -> k] becomes absolute value (data - p1)
		 * x_plus_tmp[0 -> N - k] is just data -p1 (because already sorted)
		 */

		for( int q = 0; q< N; q++){
			if (q < X_left)
				x_minus_tmp[q] = abs(data[q] - p1);
			else
				x_plus_tmp[q - X_left] = data[q] - p1;
		}

		/* start_A = 0.7
		 * q = [0 -> sizeA]
		 * inc_A = incrementAlpha
		 * alpha = 0.7 + q * inc_A
		 * Basically alpha goes from the start value through sizeA increments of inc_A 		 *
		 */


		for ( int q = 0; q < sizeA; q++){
			X_minus[q] = 0;
			X_plus[q] = 0;

			double alpha = start_A + q * inc_A;
			/* Goes through the array, If this is matlab then the code would be:
			 * alpha[sizeA] = 0.7 + sizeA*inc_A
			 * data[N] = data vector coming in
			 * k = index value			 *
			 * x_minus_tmp [k];
			 * x_plus_tmp[N-k];
			 * x_minus_tmp = abs( data[0-k] - data[k] );
			 * x_plus_tmp = data[k-N] - data[k];
			 * X_minus[q] = sum (x_minus_tmp.^alpha[q])/N;
			 * X_plus[q] = sum (x_plus_tmp.^alpha[q])/N;
			 *
			 * So the outer loop pins k, the inner loop pins q.
			 * and for each nested loop we have X_minus[k, q] and X_plus [k,q]
			 *
			 */

			for (int t = 0; t < N; t++){
				if (t < X_left)
					X_minus[q] += pow(x_minus_tmp[t], alpha) / N;
				else
					X_plus[q] += pow(x_plus_tmp[t - X_left], alpha)/ N;
			}
			/////////////////////////////////////////////////////////////////////////////////////////
			// Debugging code
			////////////////////////////////////////////////////////////////////////////////////////////
	/*		if (k == 400){
				if (q == 40){
					for (int i = 0; i<k; i++)
						cout << "host pow(x_minus_tmp["<< i << "],alpha) = "  << pow(x_minus_tmp[i], alpha) << endl;
				}
			}

			if (k == 275)
				cout << "host X_minus[" << q << "] = " << X_minus[q] << endl; */
		}

		//Produce potential parameter sets as a function of  Alpha, Sigma, Kappa which will be used for maximization
		for (int q = 0; q < sizeA; q++){

			double alpha = start_A + q * inc_A;
			Kappa[q][k] = pow( X_minus[q] / X_plus[q], 1 / (2*(alpha +1)));

			Sigma[q][k] = pow((alpha*  pow( X_minus[q]*X_plus[q], alpha/ (2*(alpha+1)))\
				* ( pow(X_plus[q], 1/(alpha+1)) + pow(X_minus[q], 1/(alpha+1)))), 1/alpha);

			double sgngam = 0;
			H[q][k] = log(Sigma[q][k]) + locallngamma(1/alpha, &sgngam) \
				+ log(Kappa[q][k] + 1/Kappa[q][k]) + 1/alpha - log(alpha);

			//cout << "Kappa : " << Kappa[q* N + k] << "\t" << "Sigma : " << Sigma[q* N + k] << "\t"<< "H :" << H[q* N + k] << endl;
			// keep record of the minimum entropy index
			if(q < sizeA -2 && q >=1 && H[q][k] < min_H ){
				min_H = H[q][k];
				idx_k = k;
				idx_q = q;
			}
		}

		delete [] x_minus_tmp;
		delete [] x_plus_tmp;
	}

	///////////////////////////////////////////////////////////////////////////////////////////
	// Debug code to print out and compare
	//////////////////////////////////////////////////////////////////////////////////////////
/*	for (int q = 0; q < sizeA; q++){
		double minH2 = 1e10;
		double tempKappa;
		double tempSigma;
		for (int k = 0; k < N; k++){
			if (H[q][k] < minH2){
				minH2 = H[q][k];
				tempKappa = Kappa[q][k];
				tempSigma = Sigma[q][k];
			}
		}
		cout << "Kappa[" << q << "] = :" << tempKappa;
		cout << "      Sigma[" << q << "] = :" << tempSigma;
		cout << endl;
	}
*/
	///////////////////////////////////////////////////////////////////////////////////////////
	// end debug
	//////////////////////////////////////////////////////////////////////////////////////////

	params[0] = data[idx_k]; // theta
	params[1] = Sigma[idx_q][idx_k]; // sigma
	params[2] = Kappa[idx_q][idx_k]; // kappa
	params[3] = start_A + idx_q * inc_A; // alpha

	for (int i = 0 ; i< sizeA ; i++){
		delete [] Kappa[i];
		delete [] Sigma[i];
		delete [] H[i];

	}

	delete [] Kappa;
	delete [] Sigma;
	delete [] H;
	delete [] X_minus;
	delete [] X_plus;
	//cout << "Got here with no problem "<< endl;
}


/////////////////////////////////////////////////////////////////////////////////////////////////
// The following code creates a matrix that does the calculations of  ( p1-data[k] ) ^ alpha
// The code creates an [i][j] = [600][80] matrix where entry D[i][j] = Div[i] ^ alpha[j]
// Div[i] = minDIV * i for i from 0 to 599
// alpha[i] = start_A + inc_A * i = 0.7 + 0.1 * i for i = 0 to 79.
////////////////////////////////////////////////////////////////////////////////////////////////

void MakeDeviationMatrix(double **D){

	double inc_A = 0.1;//////Iteration Increment of Alpha Parameter
	double start_A = 0.7;//////// Value to Start Alpha Iteration
	double alpha[sizeAlpha];
	double DIV[numDIV];
	////double **D = new double*[sizeALPHA];
	for (int i = 0 ; i< numDIV ; i++)
	{
		D[i] = new double[sizeAlpha];
	}

	for(int i = 0; i<sizeAlpha; i++)
	{
		alpha[i] = start_A+inc_A*i;
	}

	for(int i = 0; i<numDIV; i++)
	{
		DIV[i] = minDIV*i;
	}

	for(int j = 0; j<numDIV; j++)
	{
		for(int i = 0; i<sizeAlpha; i++)
		{

			D[j][i] = pow(DIV[j],alpha[i]);
		///	fprintf(dMat,"%lf\t",D[j][i]);
		}
		///fprintf(dMat,"\n");
	}
	///fclose(dMat);
}

/*void fast_AEPD_Fit(double * x, int N, double * params, struct DeviationStruct *D){

	vector<double> data(x, x+N);
	sort(data.begin(),data.end());



	double inc_A = INCREMENT_ALPHA;//////Iteration Increment of Alpha Parameter
	double start_A = 0.6;//////// Value to Start Alpha Iteration
	int sizeA = sizeALPHA;//// Number of iterations to perform

	double Kappa[sizeALPHA][NUM_FITS*WINDOW_SIZE];
	double Sigma[sizeALPHA][NUM_FITS*WINDOW_SIZE];
	double  H[sizeALPHA][NUM_FITS*WINDOW_SIZE];
	////for (int i = 0 ; i< sizeA ; i++){
		///Kappa[i] = new double[N];
		///Sigma[i] = new double[N];
		///H[i] = new double[N];
	///}

	double  X_minus[sizeALPHA];
	double  X_plus[sizeALPHA];


	int idx_k = 0;
	int idx_q = 0;
	double min_H = 10000000000;

	int stepX = round(N/NUM_STEPS);
	int X_right = 0;
	int X_left = 0;
	if (stepX <=1)
		stepX=1;
	int x_minus_tmp[NUM_FITS*WINDOW_SIZE];
	int x_plus_tmp[NUM_FITS*WINDOW_SIZE];

	double p1;

	for(int k  = EDGE -1 ; k<= N- EDGE -1 ; k=k+ stepX){

		p1 = data[k];

		X_right = N- k-1;
		X_left = k+1 ;


		for( int q = 0; q< N; q++){
			if(q<=k){
				x_minus_tmp[q] = round(FACTOR*abs(data[q] - p1));
				x_plus_tmp[q]=0;
			}
			else if(q>k){
				x_plus_tmp[q] = round(FACTOR*abs(data[q] - p1));
				x_minus_tmp[q]=0;
			}
		}


		for ( int q = 0; q < sizeA; q++){
			X_minus[q] = 0;
			X_plus[q] = 0;

			for( int t = 0; t< N; t++){
				if(x_minus_tmp[t]>0){
					///cout<<"x_minus:\t"<< D[x_minus_tmp[t]][q]<<endl;
					X_minus[q] += D[x_minus_tmp[t]].shapeIDX[q] / N;
				}
				if(x_plus_tmp[t]>0){

					X_plus[q] += D[x_plus_tmp[t]].shapeIDX[q]/ N;
				}
			}
		}


		//Produce potential parameter sets as a function of  Alpha, Sigma, Kappa which will be used for maximization
		for (int q = 0; q < sizeA; q++){

			double alpha = start_A + q * inc_A;
			Kappa[q][k] = pow( X_minus[q] / X_plus[q], 1 / (2*(alpha +1)));

			Sigma[q][k] = pow((alpha*  pow( X_minus[q]*X_plus[q], alpha/ (2*(alpha+1)))\
				* ( pow(X_plus[q], 1/(alpha+1)) + pow(X_minus[q], 1/(alpha+1)))), 1/alpha);

			double sgngam = 0;
			H[q][k] = log(Sigma[q][k]) + lngamma(1/alpha,sgngam) \
				+ log(Kappa[q][k] + 1/Kappa[q][k]) + 1/alpha - log(alpha);

			//cout << "Kappa : " << Kappa[q* N + k] << "\t" << "Sigma : " << Sigma[q* N + k] << "\t"<< "H :" << H[q* N + k] << endl;
			// keep record of the minimum entropy index
			if(q < sizeA -1 && q >=1 && H[q][k] < min_H ){
				min_H = H[q][k];
				idx_k = k;
				idx_q = q;
			}
		}

		///delete [] x_minus_tmp;
		///delete [] x_plus_tmp;
	}

	params[0] = data[idx_k]; // theta
	params[1] = Sigma[idx_q][idx_k]; // sigma
	params[2] = Kappa[idx_q][idx_k]; // kappa
	params[3] = start_A + idx_q * inc_A; // alpha

	for (int i = 0 ; i< sizeA ; i++){
		delete [] Kappa[i];
		delete [] Sigma[i];
		delete [] H[i];

	}

	delete [] Kappa;
	delete [] Sigma;
	delete [] H;
	delete [] X_minus;
	delete [] X_plus;

}*/
